function Ef = FindFermiParab(Dv,Dc,Ev,Ec,x,Vx,BackGround,kT)

% An energy range in which to find the fermi level
Emin = -5;
Emax = +5;

% Initial setup
Ef1 = Emin; Ef2 = Emax;
Ef = 0;

for i = 1:30
    % Computing the charge in real space grid
    Sx = Dv*kT*log(1+exp((Ev-Ef-Vx)/kT))-Dc*kT*log(1+exp(-(Ec-Ef-Vx)/kT));
    TotalCharge = sum(BackGround + Sx);
    % Binary Searching
        if TotalCharge > 0
            Ef1 = Ef;
            Ef = (Ef + Ef2)/2;
        else
            Ef2 = Ef;
            Ef = (Ef1 + Ef)/2;
        end
end


end